Metropolis Monte Carlo algorithm based on reparametrization invariance 
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Abstract 

We introduce a modification of tfie well-known Metropolis importance sampling algorithm by using a methodology inspired 
on the consideration of the reparametrization invariance of the microcanonical ensemble. The most important feature of the 
present proposal is the possibility of performing a suitable description of microcanonical thermodynamic states during the 
first-order phase transitions by using this local Monte Carlo algorithm. 

PACS numbers: 05.70.-a; 05.20.Gg 
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I. INTRODUCTION 

The basic problem in equilibrium statistical mechan- 
ics is to compute the phase space average, in which 
Monte Carlo method plays a very important role P, H, [1] . 
Among all admissible statistical ensembles used with the 
above purpose, the microcanonical ensemble provides the 
most complete characterization of a given system in ther- 
modynamic equilibrium. While microcanonical calcu- 
lations based on microscopic dynamics can suffer from 
the existence of metastable states with large relaxation 
times [1, [1] , the available Monte Carlo methods based on 
reweighting technics (such as the Multicanonical Monte 
Carlo [1, Q) are also so much machine-time consuming 
due to they try to obtain the density of states for all 
values of energy in only one run. It is very well known 
that the microcanonical observables can be suitably ob- 
tained by using the canonical ensemble whenever ensem- 
ble equivalence takes place in a sufficient large system. 

One of the most famous and widely used Monte Carlo 
algorithm considering the Gibbs canonical ensemble is 
the Metropolis importance sampling algorithm This 
method has some important features. It is extremely gen- 
eral and each moves involves O (1) operations. However, 
its dynamics suffers from critical slowing down, that is, 
if N is the system size, the correlation time r diverges 
as a critical temperature is approached: the divergence 
follows a power law behavior t oc in second-order 
phase transitions, while r diverges exponentially with N 
in first-order phase transitions as consequence of the en- 
semble inequivalence (multimodal character of the energy 
distribution function) . 

The aim of the present Letter is to develop a new 
methodology in which the local algorithm Metropolis 
Monte Carlo could be used for performing microcanonical 
calculations by avoiding the anomalous behaviors associ- 
ated with the presence of a first order phase transition 
in a given energetic region. The key of our method is 
to consider certain generalized canonical-like ensemble 
which is equivalent to the microcanonical ensemble for 
N large enough in spite of the presence of a first-order 
phase transition, a methodology inspired by our recent 
paper |9|. 



II. THE METHOD 

Let us consider a given Hamiltonian system H becom- 
ing extensive in the thermodynamic limit N —f oo. As 
already commented, the basic idea of our method is to 
substitute the Gibbs canonical ensemble by the following 
generalized canonical ensemble: 



1 



exp 



{-vQ ( 



H 



(1) 



where Q (E) = Nip (E/N), being ip (e) certain analytical 
bijective function of the energy per particle of the system 
e = E/N. 

The Planck potential V {rj, N) ~ —\nZc {rj, N) associ- 
ated with the partition function: 



exp 
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can be estimated for N large by using the steepest de- 
scend method as follows: 
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where s (e) = In {W {E, N)) /N is the entropy per 
particle of the system, being W — H. {E, N) 6eo the 
microcanonical accessible volume, and Se, the only 
one stationary energy of the maximization problem 
max^^ {rjif (s) — s (e)} whose stationary conditions de- 
mand: 
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We found by combining ([4]) and ([5]) that the necessary 
and sufficient condition for the existence of a unique sta- 
tionary point of the above maximization problem for any 
r] is given by: 



dpje) d j f dip{e) Y'^ dsje) 
de de \ \ de J de 



> 0, 



(6) 



1 



which has to be appHcable for all those admissible val- 
ues of the energy per particle e. This condition can be 
rephrased by considering the inverse function s (</?) of the 
bijective function ip{e), and taking the entropy s as a 
scalar function which is now rewritten in terms of the 
variable tp: 



\ de 



> 0, (7) 



which clarifies us that the entropy must be a concave 
function in this last parametrization. 

The approximation ([3]) in the thermodynamic limit 
is just the Legendre transformation between the ther- 
modynamic potentials, P{ri,N) ~ riQ — 5 (8, TV) with 
S (O, N) = S {E, N), which takes place as a consequence 
of the equivalence between the generalized canonical en- 
semble ([T]) and the microcanonical ensemble. Thus, we 
can ensure the ensemble equivalence for all energy values 
by using an appropriate bijective function ip (e) without 
mattering about the presence of a first-order phase tran- 
sition. This remark is precisely the key of the Metropolis 
algorithm described below. 

Considering the analogy with the Gibbs canonical en- 
semble, the generalized canonical ensemble represents a 
system in which the parameter rj associated with the 
macroscopic observable Q is keeping fixed by means 
of certain external experimental arrangement. The 
state above is just a consequence of the reparametriza- 
tion invariance of the microcanonical ensemble: the 
microcanonical description does not depend on the 
reparametrization Q = Nip{E/N) used for describing 
the equilibrium thermodynamical states of the system. 
We only discuss in the present work the applicability of 
this symmetry in improving the ordinary Metropolis im- 
portance sampling algorithm based on the Gibbs canoni- 
cal ensemble (hereafter MMC), but the interested reader 
can see the recent paper [91] about the implication of the 
reparametrization invariance in the searching of a topo- 
logical classification scheme of the phase transitions. 

Generally speaking, the main difference between 
the MMC and the Metropolis algorithm based on 
the reparametrization invariance (hereafter GCMMC) 
is the using of the canonical-like weight u (E) = 
exp {—rjQ {E)). Therefore, the probability p for the ac- 
ceptance of a Metropolis move is given by: 
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in which all O (■^) contributions in the exponential argu- 
ment have been dismissed. While every energy decreas- 
ing Ae < is always accepted in the MMC, such prob- 
ability decreases in the GCMMC due to the presence of 
the quadratic term in the exponential argument of ([H]), 
although such effect decreases with the N increasing. 



Obviously, when N is large enough, the expression of 
the acceptance probability of the GCMMC differs effec- 
tively from the one used in the MMC by the consideration 

of a variable canonical parameter /3 (e; 77): 
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which fluctuates around the microcanonical value 
/3 (Ce) = ds (Se) /de derived from the condition Tak- 
ing into account the sharp Gaussian localization of the 
integral ([2]) in the thermodynamic limit, the dispersion 
SP = \J (5/3^) can be estimated in terms of the energy 
dispersion be = ^ (fe^) as follows: 
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where (• • • ) denotes the average over the ensemble ([T]). 
Considering the Gaussian estimation of the energy dis- 
persion as (Je^) — a1 {e^ /N with (ffe) given in ([5]), 
the equation (jlO[) can be rewritten in the following form: 
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This expression is a very nice result which talks about 
the limits of precision in a general calculation of the 
caloric curve by using a local algorithm based on the 
generalized canonical ensemble ([T]): while temperature 
can be fixed (bf3 = 0) wherever the entropy is a con- 
cave function, the inverse temperature dispersion bp can 
not vanish wherever the entropy be a convex function 
in terms of energy e. Supercritical slowing down as- 
sociated with the MMC method based on the Gibbs 
canonical ensemble is precisely related with the down- 
fall of the Gaussian estimation of the energy dispersion 

be = [—Nd'^s (e) /de^) ^ whenever the second deriva- 
tive of the entropy goes to zero and becomes nonneg- 
ative. It is very interesting to notice that the com- 
plementary macroscopic observables energy and temper- 
ature obey the uncertainly relation bEbj3 > 1 wher- 
ever d'^S {E, N) /dE"^ > 0, which is rather analogue to 
the Quantum Mechanics uncertainly relation bEbr ~ h. 
A better analysis allows us to obtain the inequalities 

be < (— A^^^s (e) /Qe^) ^ whenever entropy s (e) be a 

1 

concave function, while bf3 > {d'^s{e) jde"^ /N^'^ when 
entropy is convex. 

As already commented, the success of the GCMMC 
methods relies on the selection of a good bijective func- 
tion <p (e) satisfying the condition ([7]), which obviously 
can be done by following different schemes. For ex- 
ample, we can rephrase the stationary condition ([4]) as 
T] = ds{ip) /dip and substitute it in The resulting 
equation expresses the monotonic decreasing of the de- 
pendence r] versus the energy per particle e: 
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which talks about a bijective correspondence between 
ry and e as a consequence of the ensemble equivalence. 
This last demand suggests us the implementation of an 
indirect Monte Carlo method based on the generalized 
canonical ensemble ([ij where an arbitrary dependence 
77 (e) satisfying the above condition be assumed a priory, 
but the unknown associated bijective function tp (e) must 
be reconstructed by using reweighting schemes analogue 
to the ones used in the Multicanonical calculations [i, [J ■ 
Nevertheless, we are interested in the present work in 
the implementation of a variant of the Metropolis algo- 
rithm with acceptance probability ([H]) where the bijective 
function if (e) will be proposed a priory in order to avoids 
the ensemble inequivalence in a Hamiltonian system with 
short-range interactions. 

Let us assume that the interest system exhibits an en- 
semble inequivalence inside the energy interval (£i,£2), 
where (ci) /9e^ < 0. The microcanonical parameter 
P (e) = ds (e) /de will change slowly with the energy e in 
this region due to |9^s(e) /de^\ ~ in the first-order 
phase transitions (see in the application example de- 
scribed in the next section), so that, 77 (e) dip (e) /de ~ Pc 
where [3^ is the inverse critical temperature. Assuming 
this last relation as a first approximation, its substitution 
in the condition of the ensemble equivalence ([5]) yields: 
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with A (e) = —9 In 77 (e) / de. This demands is very easy to 
satisfy by considering the function A (e) as a large enough 
positive constant, A (e) = A > 0. This assumption allows 
us to propose the first derivative of the bijective function 
(fi {e) as follows: 



He) 
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1 if (e) > £2, 

exp (-A (£2 - e)) if (e) e (ei, £2) , 
exp (—A (£2 — £1)) otherwise, 

which should avoid the ensemble inequivalence within the 
energetic range (£i,£2). Since ^ (e) is a constant when 
(£) > £2 or (£) < £1, the GCMMC becomes in the MMC 
algorithm outside the interval (£i,£2). 

Although the piecewise function (e) is referred in 
terms of the instantaneous values of the energy e in the 
Metropolis dynamics, we propose in that the selec- 
tion of the interval during the whole computation of the 
microcanonical averages at a given energy depends on 
the average energy (£) instead of the instantaneous val- 
ues £. It can be numerically checked that the existence 
of any abrupt change in the first derivatives of the func- 
tion ^ (£) in the acceptance probability ([8]) can provoke 
a significant perturbation in the convergence of the sec- 
ond derivative of the entropy (associated with the heat 
capacity) by using the GCMMC when the instantaneous 
energy e is close to £1 or £2. However, the using of (£) 
instead of the instantaneous value e allows the function 
^ (£) to exhibit a smooth behavior beyond the interest 
energy interval (£i,£2) during the Metropolis dynamics. 



In order to obtain a caloric curve f3 versus e with ap- 
proximately M points with (£) more or less uniformly 
distributed in the interval (£i,£2), the canonical param- 
eter T] can be increased as follows: 

Vi +c if (£) > £2, 

Vr+1 = { ?7»exp(ji.A(£2 -£2)) if (£) e (£i,£2) , 

r/i + cexp (A (£2 — £1)) otherwise, 

(15) 

where c > in the /3 step outside the interest energy 
interval (£1, £2). 

The reader may notice by reexamining the equation 
([13]) that a very large value of the numeric constant A 
leads to a very small energy dispersion Se, provoking in 
this way an increasing of the inverse temperature disper- 
sion SP as a consequence of the uncertainly relation (jlip . 
An appropriate prescription of A in order to minimize 
the dispersion = ((5/3^ -I- fe^) in the caloric curve 

(£ and P in dimensionless units) is given by: 



Pc\ ^ 1. 



(16) 



Taking into consideration all the above exposed, micro- 
canonical caloric curve P (e) = ds (e) jde can be obtained 
by using our GCMMC methods as follows: 



/3 [(£) ± fc] = ,7 (£))±<5/3, 



(17) 



where 



5/3^ = ^^77^((e^(£))-(e(£))^), (18) 

(19) 



while the second derivative of the entropy or curvature 
K (£) = d^s {e) jde^ can be given by: 



k((£}) ~ (o-eCT/J - l)/crf. 



(20) 



in accordance with the equation ()lip . The error in the 
curvature after n Metropolis iterations can be estimated 
by the formula: 



5K~(2|K|-Hl/a2)W— , 



(21) 



where r is the decorrelation time, that is, the minimum 
of Monte Carlo iterations necessary to generate effec- 
tively independent, identically distributed samples in the 
Metropolis dynamics. 



III. A SIMPLE APPLICATION 

Let us apply the GCMMC method in order to perform 
the microcanonical description of a Potts model [l3| : 
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FIG. 2: Comparative study between the GCMMC and the 
Swendsen-Wang (SW) cluster algorithm. The SW algorithm 
is unable to describe the microcanonical states with a negative 
heat capacity. Notice the bimodal character of the energy 
distribution function at /? = 1.42. 



FIG. 1: Microcanonical caloric curve /3 (e) — ds (e) /de and 
curvature k (e) = d'^s (e) /9e^ of the 25 x 25 Potts model with 
q = 10 states with periodic boundary conditions obtained by 
using the GCMMC algorithm. Notice the energetic region 
with a negative heat capacity. Errors are smaller than the 
symbols linear dimension. 



on a two dimensional lattice (here with periodic bound- 
ary conditions) of = L x L spins with q ~ 10 
possible values (components). The sum is over pairs 
of nearest neighbor lattice points only and Ci is the 
spin state at the i-th lattice point. Generally speak- 
ing, this model system admits a ferromagnetic interpre- 
tation by introducing the bidimensional vector variables 
Si = [cos (KqCTi) , sin (KqCTi)] with Kq — 2T:/q, and defining 
the total magnetization as follows M = ^i- 

As elsewhere shown, the g = 10 states Potts model ex- 
hibits a first-order phase transition associated to the en- 
semble inequivalence, which provokes the existence of a 
supercritical slowing down during the ordinary Metropo- 
lis dynamics based on the consideration of the Gibbs 
canonical ensemble. Clusters algorithms, like Swendsen- 
Wang or Wolf algorithms [ll|, do not help in this case 
[T^ due to they are still based on the consideration of the 
canonical ensemble, and consequently, the supercritical 
slowing down associated to the ensemble inequivalence 
persists [isj . 

This difficulty is successfully overcome by using the 
Multicanonical Monte Carlo method Q, which re- 
duces the exponential divergence of the correlation times 
with respect to system size to a power at the first-order 



phase transitions [13] ■ Multicanonical ensemble flattens 
out the energy distribution, which allows the computa- 
tion of the density of states O {E, N) for all values of 
E in only one run. This feature of the Multicanonical 
ensemble becomes a disadvantage when we are also in- 
terested in the direct computation of the microcanonical 
average of other microscopic observables at a given value 
of the energy, i.e. the magnetization density dependence 
m (e) = (M) /N of the Potts model ^ with q = 10. 
This aim can be easily performed by using the GCMMC 
algorithm. 

A preliminary calculation by using the MMC allows 
us to set the interest energetic window with ei — 0.2 
and £2 — 1-2 in which takes place the critical slowing 
down associated with the existence of a first-order phase 
transition in this model for L = 25. The inverse critical 
temperature was estimated as /3c — 1.4, allowing us to 
set A — 0.8. Besides, we also set M = 50 and c — 0.02 in 
the rj increasing given by the equation (|15p . The caloric 
/3 (e) and curvature k (e) curves obtained by using the 
GCMMC algorithm with n — 10^ Metropolis iterations 
for each point is shown in the FIGdJ A comparative 
study between the present method and the Swendsen- 
Wang clusters algorithm [l^ is shown in the FIGlH 

The backbending in the caloric curve (3 versus e is 
directly associated to the existence of a negative heat 
capacity c (e) — — (e) /k (e) when e E {ea,Sb) with 
Ea — 0.51 and — 0.93, which is a feature of a first- 
order phase transition in a small system becoming ex- 
tensive in the thermodynamic limit. Since the heat ca- 
pacity is always positive within the canonical ensemble. 
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FIG. 3: Magnetic properties of the model within the micro- 
canonical ensemble in which is shown the existence of two 
phase transitions at £// ~ 0.7 (ferro-ferro) and e/p — 0.8 
(ferro-para) . 
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FIG. 4: Histogram of the projection of the microscopic mag- 
netization along the direction of the spontaneous magnetiza- 
tion of the ferromanet type I, which demostrates the existence 
of metastable states when e — 0.844. 



Cc = (Ae^)^ > 0, such anomalous regions are in- 
accessible in this description {{■ • denotes the canoni- 
cal average). This fact evidences the existence of a sig- 
nificant lost of information about the thermodynamical 
features of the system during the occurrence of a first- 
order phase transitions when the canonical ensemble is 
used instead of the microcanonical one. This difficulty is 
successfully overcome by the GCMMC algorithm, which 
is able to predict the microcanonical average of the mi- 
croscopic observables in these anomalous regions where 
any others Monte Carlo methods based on the considera- 
tion of the Gibbs canonical ensemble such as the original 
MMC, the Swendsen-Wang and the Wolff single cluster 
algorithms certainly do not work. 

This fact is clearly illustrated in the FlGl2l which 
evidences that the Swendsen-Wang algorithm is unable 
to describe the thermodynamic states with a negative 
heat capacity: its results within the anomalous region 
differ significantly from the ones obtained by using the 
GCMMC algorithm. The Swendsen-Wang dynamics ex- 
hibits here an erratic behavior originated from the com- 
petition of the two metastable states present in the 
neighborhood of the critical point (the energy distribu- 
tion function in the canonical ensemble is bimodal when 
P e (/3i,,32), where Pi = 1.405 and /?2 = 1.445, being 
this feature the origin of the supercritical slowing down) . 

As already shown by Gross in ref. 15], a negative heat 
capacity in a system with short-range interactions can be 
associated to the existence of a non- vanishing interphase 
surface tension. Other quantities like the transition tem- 
perature (3cr and the latent heat qiat can be derived easily 
from the caloric curve I3{e). Following the same proce- 



dure developed by Gross, our calculations allow us to 
obtain the following estimates for the g = 10 states Potts 
model with L = 25: (3cr ^ 1-421 and qiat 0.78. 

The study of the magnetic properties of this model 
within the microcanonical ensemble illustrated in the 
FIG131 shows us that the existence of the non-vanishing 
interphase surface tension is not the only one informa- 
tion hidden behind of a negative heat capacity within 
the canonical description. The modulus of the magneti- 
zation density m (e) evidences what could be considered 
as the signature of two phase transitions within the mi- 
crocanonical description of this model system: a continu- 
ous (paramagnetic-ferromagnetic) phase transition at the 
critical point £fp — 0.8, and a discontinuous (ferro-ferro) 
phase transition at £// ~ 0.7 (from a low magnetized 
ferromagnetic type II phase towards a high magnetized 
ferromagnetic type I one). The investigation of these in- 
teresting behaviors deserves a further study. 

Most of thermodynamic points of this dependences 
were obtained from a data of n = 10^ Metropolis it- 
erations, with the exception of all those points belong- 
ing to the energetic interval (0.7,0.93) where very large 
fluctuations of the magnetization density were observed. 
The possibility of the GCMMC algorithm of performing 
a localized computation of the microcanonical averages 
at a given value of the energy allows us to increase the 
Metropolis iterations n up to 5 x 10^ for each point within 
the above energetic range in order to reduce the signifi- 
cant dispersion of the expectation values observed there. 

The very large fluctuations, the peculiar behavior of 
the dispersion g (e) = <^(M- (M))^^ /N, the qualitative 
form of the magnetization curve m (e) shown in FIG131 
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as well as the apparent large relaxation times of the ex- 
pectation values during the GCMMC dynamics, suggest 
us the presence of several metastable states with different 
magnetization densities at a given energy within this last 
region, whose existence is shown in the FIG HI 

The anomalies observed in the thermodynamic char- 
acterization of the q = 10 state Potts model suggest 
strongly the existence of the critical slowing down phe- 
nomena within the " microcanonical description" pro- 
vided by the GCMMC algorithm. This fact inspires the 
development of nonlocal Monte Carlo algorithms based 
on the consideration of the generalized canonical ensem- 
ble ([1]) in order to address the microcanonical description 
of more larger systems undergoing a discontinuous (first- 
order) phase transitions within the canonical description. 

IV. CONCLUSIONS 

We have presented a very simple methodology for im- 
proving the ordinary Metropolis importance sampling al- 
gorithm 8j in order to allow the computation of the mi- 
crocanonical averages during the occurrence of the first- 
order phase transitions in systems with short-range in- 



teractions. The key of our approach is the consideration 
of certain generalized canonical ensemble ^ inspired on 
the reparametrization invariance of the microcanonical 
description [9^, which becomes equivalent to the micro- 
canonical ensemble when the size N of the interest system 
is large enough. 

A direct advantage of this new methodology is the pos- 
sibility of performing a localized computation of the any 
microcanonical average at a given value of the energy, 
and avoid in this way the unnecessary computation of 
the whole energy range, a feature of many other Monte 
Carlo technics based on reweighting the energy histogram 
[as , [l^ . The existence of anomalies within the micro- 
canonical description of the q = 10 states Foots model 
inspires the further development of nonlocal Monte Carlo 
algorithms based on the consideration of the generalized 
canonical ensemble ([l}. 
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